home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 209_01 / simplexr.doc < prev    next >
Text File  |  1990-03-03  |  10KB  |  335 lines
  1. SIMPLEXR.DOC       VERS:- 01.00  DATE:- 09/26/86  TIME:- 10:02:52 PM
  2.  
  3. notes on the simplex method of function minimization, 
  4.     by use of the Nelder-Mead algorithm.
  5.  
  6.  
  7.  
  8.  
  9.  
  10.  
  11.         By J. A. Rupley, Tucson, Arizona
  12.  
  13.  
  14.          NOTES ON THE ALGORITHM OF THE SIMPLEX FITTING PROGRAM
  15.  
  16.  
  17.         A. THE SIMPLEX METHOD OF FUNCTION MINIMIZATION
  18.  
  19.  
  20.              The fitting of a model with variable parameters to a set of
  21.         data points is handled by minimizing an error function, such as
  22.         the sum of squares of the differences between observed values and
  23.         values calculated according to the model.
  24.  
  25.              We describe here a simplex approximation procedure for
  26.         estimating the parameter values that give a function minimum.
  27.         The strategy is that of Nelder and Mead, and their article
  28.         (Computer Journal, vol 7, p 308, 1965) should be consulted for
  29.         more information than is given below.
  30.  
  31.              A simplex is constructed in parameter space, with vertices
  32.         that are arbitrary but not too far from the point at which the
  33.         function is a minimum and that describe a volume that includes
  34.         the minimum point.
  35.  
  36.              The position of each vertex in parameter space defines (or
  37.         equally, is defined by) a set of parameter values, for which one
  38.         calculates the function value for that vertex.  The simplex has
  39.         (M + 1) vertices, where M is the number of free parameters.
  40.  
  41.              In successive iterations the vertex with the highest
  42.         function value is moved, to obtain a new vertex position of
  43.         smaller function value.  The movement is directed with respect to
  44.         the center of the reduced simplex, which is the simplex less the
  45.         highest vertex.  The movements can be the following:
  46.              reflection;
  47.              expansion beyond the reflected position;
  48.              reflection after failed expansion;
  49.              contraction toward the center of the reduced simplex from
  50.                   the original position;
  51.              contraction from the reflected position;
  52.              if none of these operations gives a lower function value
  53.                   then the entire simplex is shrunk about the lowest
  54.                   vertex.
  55.  
  56.              Exit from the iterative minimization is if the fractional
  57.         RMS deviation of the function values at the vertices is less than
  58.         a test value and if the centroid of the simplex is within two RMS
  59.         deviations of the mean of the vertices.  The default test value
  60.         is 1E-8 * the mean function value.
  61.  
  62.              Thus at exit from the minimization, the centroid of the
  63.         simplex gives an arbitrarily close approximation of the parameter
  64.         values at the function minimum.
  65.  
  66.  
  67.                                         1
  68.  
  69.  
  70.  
  71.  
  72.  
  73.  
  74.  
  75.  
  76.  
  77.  
  78.              The routine that calculates the function values must be
  79.         adapted to the model to be optimized.
  80.  
  81.              Entry of the data and starting parameter values can be at
  82.         compilation, by initializing the data and parameter arrays, or by
  83.         reading a text file.
  84.  
  85.              Results should be written to a file, for convenience of
  86.         review and for recovery from crashes.
  87.  
  88.  
  89.  
  90.  
  91.  
  92.  
  93.  
  94.  
  95.  
  96.  
  97.  
  98.  
  99.  
  100.  
  101.  
  102.  
  103.  
  104.  
  105.  
  106.  
  107.  
  108.  
  109.  
  110.  
  111.  
  112.  
  113.  
  114.  
  115.  
  116.  
  117.  
  118.  
  119.  
  120.  
  121.  
  122.  
  123.  
  124.  
  125.  
  126.  
  127.  
  128.  
  129.  
  130.  
  131.  
  132.  
  133.                                         2
  134.  
  135.  
  136.  
  137.  
  138.  
  139.  
  140.  
  141.  
  142.  
  143.  
  144.         B. THE QUADRATIC FIT TO THE LEAST SQUARES FUNCTION SURFACE
  145.  
  146.  
  147.              Standard deviations of the parameters are calculated by
  148.         fitting a quadratic function to the surface about the minimum in
  149.         parameter space and then using the properties of this function to
  150.         calculate the variance-covariance matrix for the parameters.  The
  151.         strategy is a modification of that of Nelder and Mead (1965).
  152.  
  153.              The contracted simplex obtained through the previous
  154.         minimization process is used to define a new coordinate system of
  155.         M oblique axes in parameter space.  The origin of the axes is at
  156.         the centroid of the simplex.  Each axis corresponds to a free
  157.         parameter.  The unit value (scale) for each axis is set in the
  158.         new system at the average of the deviations of the vertices of
  159.         the simplex from the centroid value.
  160.  
  161.              In effect, one constructs in the new coordinate system a new
  162.         simplex, based on the one obtained by the previous minimization
  163.         process.  The (M + 1) vertices of the new construct are at the
  164.         centroid (the zero of the new system) and at the unit points on
  165.         each free parameter axis.  As before, the values of the least
  166.         squares function at these vertices describe a surface.  It is
  167.         this least squares surface to which a quadratic function is to be
  168.         fit.
  169.  
  170.              If necessary, the scales of the axes (the unit values) in
  171.         the new coordinate system can be adjusted to give better behavior
  172.         of the surface in the vicinity of the minimum.
  173.  
  174.              The diagonal matrix Q transforms the new coordinate system
  175.         back to the original system.  The elements of Q are for each
  176.         parameter the (adjusted) average differences between the vertices
  177.         and the centroid of the original simplex.
  178.  
  179.              The function value, y, in the region near the function
  180.         minimum generally can be approximated in the new coordinate
  181.         system by the quadratic vector function:
  182.  
  183.               y  =  a(o)   +   2 a'.x   +   x'.B.x                   (1)
  184.  
  185.              The elements of the vector x are the values of the M free
  186.         parameters at a point in parameter space.  The scalar a(o), the
  187.         vector a and the matrix B are determined by the shape of the
  188.         surface y near the minimum and can be calculated with the use of
  189.         simple numerical approximations, evaluated at x = 0 (the position
  190.         of the centroid of the original simplex):
  191.  
  192.  
  193.  
  194.  
  195.  
  196.  
  197.  
  198.  
  199.                                         3
  200.  
  201.  
  202.  
  203.  
  204.  
  205.  
  206.  
  207.  
  208.  
  209.                 a(0)   =   y(0)
  210.  
  211.                 ai     =   0.5 * (dy/dxi)
  212.                        =   0.25 * (yi - y-i)
  213.  
  214.                 bij    =   0.5 * (d2y/dxi.dxj)
  215.                        =   0.25 * (yij + y-i-j + 2 * y(0)
  216.                                            - y-i - y-j - yi - yj)
  217.  
  218.                 bii    =   0.5 * (d2y/dxi2)
  219.                        =   0.5 * (yi + y-i - 2 * y(0))
  220.  
  221.              The point x = 0 (the centroid position) may not be the true
  222.         position of the function minimum, although it is assumed to be
  223.         close to it.  A refined estimate of the minimum position, based
  224.         on the quadratic approximation of equation (1) above, is given by
  225.         x(min):
  226.  
  227.              xmin = - B-1.a
  228.  
  229.  
  230.              The function minimum at x(min) is given by y(min) :
  231.  
  232.              ymin = a(0) - a'.B-1.a
  233.  
  234.  
  235.              The position of the minimum in the original coordinate
  236.         system is obtained by back transformation of x(min) with the
  237.         matrix Q, giving the point p(min):
  238.  
  239.              pmin = p(0) - Q'.B-1.a
  240.  
  241.              As a test of the quality of the quadratic approximation (1),
  242.         there should be close agreement between the function values
  243.         y(min), y(p(min)), and y(centroid), and between the sets of
  244.         parameter values defining the points p(min) and the centroid.
  245.  
  246.              The variance-covariance matrix C is given by
  247.  
  248.              C = Q'.B-1.Q
  249.  
  250.              A diagonal element of C, multiplied by the mean square error
  251.         (MSE), gives the variance of the corresponding parameter:
  252.  
  253.              var(i) = cii * MSE
  254.  
  255.              MSE = ymin/(n - m)
  256.  
  257.         where the divisor (n - m) = # degrees freedom: n = # of data
  258.         points used in calcn of ymin, and m = # of free parameters; and
  259.         ymin = the sum of residuals squared
  260.  
  261.  
  262.  
  263.  
  264.  
  265.                                         4
  266.  
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  
  273.  
  274.  
  275.              The standard deviation is the square root of the variance:
  276.  
  277.              std-dev(i) =